clear all;
%-------------------------------------------------------------
% runfem2d version 1.0
%--------------------------------------------------------------
% Copyright 2010 Oka Kurniawan (kurniawano@ihpc.a-star.edu.sg)
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%-------------------------------------------------------------

%debug_on_error(1);
%read elements of the domain
[n1,n2,n3]=textread('ushapesharp.msh.7','%d %d %d');

%read top surface boundary 
topnodes=textread('ushapesharp.msh.8.p','%d');

%read surface boundary for pn side surface
pnnodes=textread('ushapesharp.msh.9.p','%d');

%read locations for nodes
[x,y,z]=textread('ushapesharp.msh.nodes','%f %f %f');

nodesxyz=[x,y];

Nelems=length(n1);
Nnodes=length(x);

%set diffusion length
difflen=2.76; 
lambda=1/difflen;

%preprocessing from FEM solver
alphax=-ones(Nelems,1);
alphay=-ones(Nelems,1);
beta=-lambda.^2*ones(Nelems,1);
f=zeros(Nelems,1);

elemnode=[n1,n2,n3,alphax,alphay,beta,f];


bvaltop=zeros(length(topnodes),1);
bvalpn=ones(length(pnnodes),1);
bval=[bvalpn;bvaltop];
bnodes=[pnnodes;topnodes];

bdrnode=[bnodes,bval];
[result,K,B,idxComp,idxBn,Kbc,retBnval]=solvefem2dvec(nodesxyz,elemnode,bdrnode);

%obtaining solutions
solution=[result;retBnval];
xs=[x(idxComp);x(idxBn)];
ys=[y(idxComp);y(idxBn)];


%xuni=(0:0.2:10)';
%yuni=(0:0.2:5)';
%[xi,yi]=meshgrid(xuni,yuni);

%zi=griddata(xs,ys,solution,xi,yi);

%testx=xi(5,:);
%testz=zi(5,:);

%logcur=log(testz);
%plot(testx,logcur);
%[p,s]=polyfit(testx,logcur,1);

%mesh(xi,yi,zi);

xvtk=[x,y];
econ=[n1,n2,n3];
data(idxComp)=result;
data(idxBn)=retBnval;
twod_to_vtk(xvtk,econ,data,'result.vtk','charge Collection Probability');


